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Abstract. A newly developed 3-D Monte Carlo model is used, in conjunction with a multi-line non-LTE radiative 
transfer model, to determine the mass-loss rate of the Wolf-Rayet (W-R) star in the massive binary V444 Cyg 
(WN5-I-06). This independent estimate of mass-loss rate is attained by fitting the observed He I 5876 A and 
Hell 5412 A line profiles, and the continuum light curves of three Stokes parameters (7, Q, U) in the V band 
simultaneously. The high accuracy of our determination arises from the use of many observational constraints, and 
the sensitivity of the continuum polarization to the mass-loss rate. Our best fit model suggests that the mass-loss 
rate of the system is Mwr = 0.6 (±0.2) x 10~^ Mq yr~^, and is independent of the assumed distance to V444 Cyg. 
The fits did not allow a unique value for the radius of the W-R star to be derived. The range of the volume filling 
factor for the W-R star atmosphere is estimated to be in the range of 0.050 (for _Rwr = 5.0 -Rq) to 0.075 (for 
RwB. ~ 2.5Rq). We also found that the blue-side of Hel 5876 A and Hell 5412 A lines at phase 0.8 is relatively 
unaffected by the emission from the wind-wind interaction zone and the absorption by the O-star atmosphere; 
hence, the profiles at this phase are suitable for spectral line fittings using a spherical radiative transfer model. 
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1. Introduction 



The short-period {P = 4.212 days, [Khaliullin et al. 1984]) 
eclipsi ng, massive bhiaiy V444 Cyg (WN5+O0 III-V) lias 



subject of exLensive sLudies since its discuveiy 



been t 

(e.g., I ^^ilsuu 1939 | , pVlLiuch 1950 | , I g faonr- 



there are at least 6 different methods of determining the 
mass-loss rate of a W-R star: 1. dynamical method using 
the change in the orbital period of binary (e.g., KhaliuUin 
et al. 1984 ), 2. polarization variation method (St-Louis 
3. X-ray spectra synthesis method using a 



et aL 19881 ), , ..^ 

. . muwii ±aoo| , hydrodynamical model (e.g., fetevens et al. 1996; Pittard 

Marchojiiku el al. 1994rt§§7). V444 Cyg exhibits vailabil- — 5 — nnnnl a j- /tS 1^ a — — j 

t 1' I & Corcoran 2002), 4. radio/lR continuum nux method 

ity in pulaiizaliuii, liiie sLieiiglli and X-iay flux as a fuiic 



tion of orbital phase. The variability arises from occulta- 
tion of the photosphere, from perturbations induced in the 
extended atmosphere of the Wolf-Rayet star (W-R) by the 
O star and its wind, and from the wind-wind interaction 
region. Despite the complexities, many authors ([Hamann 



fe Schwarz 1992; St-Louis et al. 1993; Marchenko et al 



1994 



Cherepashchuk et al. 1995; Moffat & Marchenko 



1996; l ^laielieiiku el al. 1997 r SLeveiis IIuwaiLh W mj 



have u icd this objeet to detmmiiic fundameiiLal paiaiiie- 
ters of ' the W-R slai. 



(e.g. Wright fc Barlow 1975), 5. radiat ive transfer method 
(e.g., Hillier 198£; Schmutz et al. 1989|) and 6. photo metric 
variability method (e.g., Lamontagne et al. 1996). The 
mass-loss rate estimated from methods 4 and 5 are, in 
general, about 2 — 3 times higher than the values estimated 
from methods 1, 2 and 3. 

In the case of V444 Cyg, the mass-loss rate mea- 
sured from the orbital period change by KhaliuUin 

(19901) 



(19841 ), [UnderhiU et al. 



One of the most important and uncertain parameters 
in the stellar evolution calculation of a massive star is the 
mass loss rate (M) during the W-R star phase. Presently 



et al 

et al. (1995D are Mwr = 1.0 
-1 



and lAntokhin 

1 



-^Mc:, yr- 
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lO"'^ Mq yr" 

0.4 X IQ-^AfQyr-i and Mwr = 0-7 x 10 
respectively. (See also Table |l|.) [Prinja et al. (1990|) 
found A^wR — 2.4 x 10~^MQyr^^ from the free-free 
radio emission flux. By modeling the infrared lines, 

= (2.5 - 5) X 
(19981) pub- 



Howarth 



10-=' M, 



Schmutz (1992 ) obt ained A/wr - 
^. More recently, Nugis et al. 



2 
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lished the "dumping-corrected" mass-loss rates for 37 
W-R stars from the radio emission power and the spec- 
tral index {a = din f,^/dliii'). They obtained Mwr = 
0.92 X IQ-'^ M(7^ yi-^ for V444 Cyg. 

Hamann & Schwarz (iyyz'), by using a non-LTE ra- 
diative transfer model with a spherically expanding atmo- 
sphere, simultaneously fitted a set of helium emission lines 
and the light curve of V444 Cyg. The mass-loss rate esti- 
mated from their analysis is Mwr — 1-26 x 10^^ Mq yi"-^. 
Similar to Methods 4 and 5 mentioned above, the model 



The time dependent spectrum of the helium lines ob- 
tained by [Shore fc Brown (1988|), [Underhill et al. (1988D, 



Marchenko et al. (1994), Marchenko et al. (1997) and 



Stevens fc Howarth (1999| ) shows very complex variabil 
ity associated with the orbital motion of the stars and 
the bow shock. Although the terminal velocities of the 
O star (2540kms-i: [Underhill et al. 1988D and the W- 
R star (1785 kms^^: Prinja et al. 1990) are compara- 
ble, the W-R star has a much denser wind, c.f., Mq w 



assumt s no clumps in the stellar wind of the W-R star; 
hence, the mass-loss rate value is most likely overesti- 
mated. An important conclusion of their work is that the 
light curve is very sensitive to the inclination angle which 
must be treated as a fitt ing parameter. 

St-Louis et al. (1988, hereafter STLl) derived an an- 
alytical expression for the mass-loss rate of a W-R star 



0.6 X lO-'^MQyr-i ( [Leitherer 1988 
"T997[) and Mwr « (0.4 



in a binary system, based on the model of Brown et al 
1978[ ) which predicts the continuum polarization from a 



binary system as a function of orbital phase. The origi- 
nal model of Brown et al. (1978|) assumes point sources 



and a f ingle scattering atmosphere (optically thin). In the Pollock 1994 



Marchenko et al 

1.26) X 10-5 Mgyr-i (Table [l[). 

The momentum carried by the W-R wind is several times 
larger than that of the O star wind, and the bow shock 
produced from the colliding stellar winds is folded toward 
the O star. Further, the orbital speed of the binary is a sig- 
nificant fraction of the terminal velocity of the W-R stellar 
wind (V^ib/l^o ^ 1/4 ); hence, the shape of the bow shock 
is not expected to be cylindrically symmetric around the 
axis joining the center of the W-R star and the O star. In 
addition, the bow shock is not smooth, i.e., th e shape of 
the shock front is distorted by in s tabilities (see Stevens fc 



same paper, STLl estimated Mwr of V444 Cyg to be 
0.9 X 10^'^ M(i^ yr~^ assuming Vqo — 2500kms~\ Later 



Gayley et al. 1997; Pittard 1998 



Pittard fc 



Stevens 1999 ) — the Kelvin- Helm holtz and the non-linea r 



St-Louis et al. (1993, hereafter STL2) corrected this value 



to Afyv R = 0.6 X 10 ''M(^yT ^ using the new estimate (I996[) 



thin-shell instab il ities discussed in Chandrasekhar (1961 ), 
Vishniac (19"83[) , [Vishniac (1994[) and [Blondin fc Marks 



of Ko = 1785 kms~" ( |Prinja et al. 1990| ). STL2 also de- 
rived analytical expressions of Q fc U Stokes parameters 
as a function of orbital phase especially near the secondary 
eclipse where the W-R is eclipsed by the O star. The ob- 
served polarization curves were fitted with their model, 
and the mass-loss rate of the W-R star in V444 Cyg was 
estimated (0.75 x lO^^MQyr^^). Although these models 



are veiiy simple and puLcuLially powerful, Llie vahdiLy Is 
still in (|LiesLiuu since Lhe models assume LliaL Llie W-R 
star has a smooth wind and its atmosphere is optically 
thin. In addition, the derived expression of AfwR in STLl 
and Q & [/ in STL2 involves an integral which diverges. To 
avoid the infinity, they fixed the lower radial integration 
limit to ~ 2_R*. 

Another important, but still uncertain parameter for 
V444 Cyg is the luminosity ratio, q = Lwr/^o 



The first X-ray evidence of the colliding winds (CWs) 
in V444 Cyg was found by [Moffat et al. (1982| ) from the 
flux variability seen in Einstein Observatory data. In addi- 
tion, Corcoran et al. (199(:) and Maeda et al. (199"9| ) found 
a similar X-ray variability using ROSAT and ASCA re- 
spectively. Theoretical aspects of CW related X-ray vari- 
ability have been d e veloped by, fo r example, Luo et al. 
1990[), [Usov (1990|), [Usov (1992|), Sevens et al. (1992[) , 



Myasnikov fc Zhekov (1993[) and [Pittard fc Stevens (1997[) . 



Beals 



(1944 ) spectroscopically estimated q for the wavelength 
range 4000 - 5000A. He obtained q = 0.21 using the emis- 
sion lines, and q = 0.19 using the absorption lines. For 



the same system, Cherepashchuk et al. (1995) estimated 
the value of g in a similar manner as Beals (1944), but 
with some modification to the calculation of the equiva- 
lent width of the emission/absorption lines. They found 
q = 0.60 ± 0.06 for AA4000-6000; however, the values 
from different lines show a large scatter, ranging from 
0.36 to 0.80. From the two principal emission lines used in 
these studies, they found q = 0.51 using Hen 5412 A and 
q = 0.59 using Hei 5876 A (see their Table 1). Lastly, 
Hamann fc Schwarz (199^ ) estimated q — 3.1 according 
to their models of the light curve and helium spectrum, 
but the value is most likely incorrect since the effect of 
clumping was not included in their model. 



In a previous paper (Kurosawa & Hillier 2001, here- 
after. Paper I), we developed a 3-D Monte-Carlo radiative 
transfer model which predicts the level of continuum and 
line polarization produced by scattering of light in a stel- 
lar atmosphere with arbitrary geometry. The model can 
predict the variability features associated with orbital mo- 
tion of a binary system, including the polarization level, 
the flux level and line profile shapes. The model can treat 
a finite sized stellar disk, multiple scattering, line absorp- 
tion of the continuum photons and emission from multiple 
light sources. To achieve high precision with a minimum 
of data points (i.e. to save computer memory), an 8-way 
tree data structure cr eated via a "cell-splitting" method 
(e.g. Wolf et al. 1999 ) is used in the model. This method 
is essential since the model assumes no symmetry in the 
atmosphere. A logarithmic grid, which is commonly used 
in spherical and axi-symmetric codes, is not readily im- 
plementable in a model of arbitrary geometry. 

The goal of this paper is to apply the 3-D Monte-Carlo 
model developed in Paper I to V444 Cyg, and to interpret 
the observed variability features (seen in the continuum 
polarization and line profile shape) as a function of or- 
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bital phase. We estimate the mass-loss rate of the W-R 
component in the binary, the orbital inclination, and the 
luminosity ratio of the two stars, by fitting the observed 
He I 5876 A and Hen 5412 A line profiles, and the con- 
tinuum light curves of three Stokes parameters (I, Q, U) 
in the V band consistently. Unlike the radio/IR/spectral 
methods discussed earlier, the derived mass-loss rate de- 
termined by the polarization method is insensitive to the 
amount of clumps in the stellar wind. The continuum po- 
larization level is proportional to the electron number den- 
sity while the thermal emissivity (of radio/IR) is propor- 
tional to the square of the density. With a compHcated 
density distribution, a Monte Carlo simulation provides 
the only method to realistically predict variability in po- 
larization, continuum flux and line strength. An original 
aim was also to obtain improved estimates of the stellar 
parameters such as the radii of the W-R star and the O 
star, the monochromatic luminosity ratio of the two stars 
and the orbital inclination angle. Due to both model and 
observation uncertainties, this proved not to be feasible. 

In we discuss our models, examine the variability 
of He I 5876 A and Hen 5412 A lines, and determine the 
mass-loss rate of the Wolf-Rayet star in V444 Cyg. The 
discussion on the results will be given in and the con- 
clusion in 

2. Analysis 

The five most important parameters which we initially 
hoped to constrain in our analysis were: the mass-loss rate 
of the W-R star (A^vvr): the luminosity of the W-R star 
(-^wr), the radius of the W-R star {Rwr), the volume 
filling factor (/) and the monochromatic luminosity ra- 
tio of the W-R and O stars, q (A) = ^WR (A) /Lq (A), at 
A = 5630A. The volume filling factor is defined as the frac- 
tional volume which contains material (the higher density 
regions or clumps), and it controls the amoun t of clumps 
in the atmosphere (see e.g., Abbott et al. 1981 ). It can also 
be thought of as the fractional length along any random 
line of sight which contains clumps. 

Since the procedure of the model fitting is rather com- 
plicated, we summarize the steps for determining the pa- 
rameters below: 

1. Understand the basic behavior of the variability seen 
in He I 5876 A and Heii 5412 A lines to determine 
at which binary phase the spectrum is least affected 
by the bow shock and by the presence of the O star 
atmospheric absorption. In other words, determine at 
which phase the spectrum is modeled best by a spher- 
ical radiative transfer model (CMFGEN). 

2. At the phase chosen to be the best in Step 1, fit 
the observed helium spectrum for different W-R star 
radii, and find the corresponding value of q. In general. 



^WR 



2 X 10^ Lq is usedQ. 



^ This luminosity value is chosen fror a the mass-luminosit y 
relation of pure Heli um stars given in Meynet et al. (1994|) , 
and MwR ~ 10 Mq ( |Marchenko et al. 1994 ). 



3. Fit the continuum I, Q, U light curves around the sec- 
ondary eclipse {(/> « 0.5) to constrain the radius of the 
W-R star, (q and Rwr o,f£ constrained.) 

4. With RwK and q values determined in Step 3, re-fit 
the helium spectrum with different values of M . 

5. Compute the continuum /, Q, U light curves for the 
^WR, I 9 and Mwr values used in Step 4, and fit only 
the / light curve by adjusting the orbital inclination 
(i) and Rq- 

6. Now there are several models which fit the observed 
/ light curve and the helium spectrum simultaneously, 
but only a certain value of Mwr will fit the Q and 
U light curves; therefore, the value of A/vvr is con- 
strained, {i, i?o o-nd Mwr are constrained.) 

2.1. Line Variability and Geometric Configuration of 
the Binary 

Since the 3-D Monte Carlo polarization model introduced 
in Paper I is not a self-consistent model (i.e., the radiation 
does not modify the opacity of the atmosphere), the opac- 
ity and the emissivity must be supplied as the inputs of a 
model. There are three possible sources of the input field: 

1. parametrized functions from a semi-analytical model, 

2. output from a spherical non-LTE radiative transfer 
model with some modifications, and 3. output from a hy- 
drodynamical calculation in conjunction with method 2. 
Method 1 is useful for testing models, but not suited for 
modeling a real object. In our analysis, method 2 is mainly 
used. 

The first step, in method 2, is to compute the ob- 
ject's spectrum using a non-LTE radiative transfer code 
(CMFGEN) which includes thousands of lines to treat the 
line-blanketing effect (Hillier & Miller 1998). For example, 



see [Hillier fc Miller (1999D and [Herald et al. (200lD for a 
detailed discussion of W-R star spectrum models across 
a wide wavelength range. CMFGEN provides the opacity 
and the emissivity of the W-R star and the O star sepa- 
rately, and is designed for a single star with a "spherically" 
expanding atmosphere; hence, the output fields (opacities 
etc.) are functions of radius only. The second step is to 
compute the variability of the Stokes parameters (/, Q, 
U) and the line profiles by using the 3-D Monte Carlo 
model as a function of the binary phase. 

The velocity field in a CMFGEN model is assumed to 



have the following form (Hillier & Miller 1999) 



V (r) 



V0 + V1 + V2 



1 + {Vo/Vcore - 1) exp [- (r - R,) /h,s] 



(1) 



where Vi 



{Voc-V,^t-Vo)il-R,/r)''\ V2 



Vc^t — R*/r)^^ , hes is an isothermal effective scale 
height in the inner atmosphere, Voo is the terminal ve- 
locity, Vcorc is the expansion velocity of the core, and 
Vc^t & /32 are the parameters for the outer parts of the 
wind. /3i is equivalent to the (3 in the classic beta-velocity 
law, V (r) = Vgo (1 — R^, / r)^ , for the "inner part" of 



bhe stellar wind. See [Najarro et al. (1997[) and [Hillier fc 
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Table 1. Basic Parameters of V444 Cyg. 



PUBLISHED VALUES 



MwR [x 10~ MqJV' 
q h LwR (A) /Lo (A)] 

-RwR \_Rq\ 
Ro [Rq] 
i [deg.] 

Mi*' (O+W-R) 
m„ (O+W-R) 
DM 

Mass (Q:W-R) \Mq] 



o.et"), i.o^'\ o.4(=), {o.stliy 0.7W, 0.9^ 

2.4'*=', 2.5-5.0'-''', 1.26<«', 0.92<°' 

0.6 ± COee^' for AA4000 - 6000 A 

~ 0.2<'' AA4057 - 6563 A 

0.25 - 0.35(-'' for AA4400 - 10000 A 

2.9(0^ < 4("', 6.3l'»\ 3.1 -5.2(''' 

8.5 ± l'"', lO.O'''', 3.76<s', 8.4 - 9.3'"' 

82.8 ± 0.9''', 78.8 ± 0.5'™', 80.8 ± 1.6'"', 83.47'"' 

-5.7'"' , -4.2'9' 

8.27'"' 

11.01 ± 0.20'"', 11.27 - 11.47'^*', 9.3'"' 

25 : 10'^', 27.9(±3.2) : 9.3 (±1.0)'"' , 37.5 : 11.3'*' 



(a) |St-Louis et al. (19931), jb) [Khaliullin et al. (1984|), (c) |UnderhilI et al. (1990|), (d) [Marchenko et al. (1997|), (e) |Prinja 



et al. (]|990|). ( f) [Howarth fc Schmutz (1992[ ), (a) |Hamann fc Schwarz (1992| ), ( h) |Clierepaslichuk et al. (1995|). (i) [Bcals (194^ 
(j) |Kuhi7l968l). (fc) |Ctierepashcliuk ( 19751), (/) Piirola fc Linnaluoto (198§f7^[m) iRobert et al. (1990l), (n) [Lundstro in fc Stenholm 



Moffat & 



(1984| ), | (o) jNugis et al (1998|), (jp) iMiinch (1950| ), (q) [Marchenko et al. (1994|) (r) jCherepashchuk et al. (1984|), 

Marche Qko"(1996|) , (t) [Underhill et al. (1988| ),1m) [Forbes et al. (1992|) , (u) |Antokhin et al. (1995]) , (w) |Rodrigucs fc Magalhaes 
(I995Q 



(*) assumed the distance to be 1.7 kpc and Et-v = O.i 
be (b-«)„ = -0.30. 



Lundstrom & Stenholm (1984) - They assumed the intrinsic color to 



Miller (1999) for more explanations on this velocity law. 
Although the line-driving of 0-star winds is well under- 
stood, the multiple-scattering of W-R winds complicates 



theoretical models, and the standard CAK theory (Castor 



Abbot' , & Klein 1975) does not apply for W-R stars. 

Fig. |l| shows the basic configuration of the model ge- 
ometry. The W-R star with its extended atmosphere is 
located at the center. Surrounding the O star, the bow 
shock due to the CWs has a paraboloidal shape, and points 
away from the W-R star since its wind is much stronger 
than that of the O star. The head of the shock front 
is tilted from the purely radial direction because of the 
rapid orbital motion. Using the ratio of the terminal ve- 
locity of the W-R component {Voo = 1785 kms~^ 
990| ) and the orbital speed (T4rb 



et al. 



Prinja 
460 kms-^)ij, 

the tilt angle of the bow shock region can be estimated as 
5 ~ tan~^(T4i.b/14o) ~ 15°. Extensive discussions on the 



morph ology of the bow shock can be found in [Marchenko 
et al. ([993), including the issue of whether the W-R wind 
is impacting onto the O star s urface or not. (See also t he 
hydrodynamical calculation of Pittard fc Stevens 1999| ). 

At 4> (orbital phase) — 0, the W-R star is in front of 
the O star, and an observer is looking at the system from 
the top of Fig. since the inclination angle i is about 
80° (almost edge-on view). At (p = 0.5, the O star is in 
front of the W-R star. The directions of an observer at 



We simply estimated the value by assuming a circular orbit 



If the radial velocity measurements of Marchenko et al. (1994 ) 
and the orbital inclination i = 78.8° ( Robert et al. 1990| ) were 
used, we would obtain Voih ~ 440 kms"'^. Consequently, the 
tilt angle would be 5 « 14°. 



(p = 0, 0.25, 0.5, 0.75 are indicated in the same figure, 
and they infer the sense of the orbital motion (counter- 
clockwise in the figure). The strong stellar wind from the 
W-R star is interrupted by the shock front, and the density 
behind the shock is assumed to be negligibly small for 
simplicity. In other words, we do not include the O star 
wind explicitly, and we place the bow shock zone rather 
artificially. 

Next, we propose a simple model of the gas flow in the 
CW zone and the geometry of the CW zone. We assume 
the shape of the CW zone to be paraboloid^ with thickness 
d. The center of the paraboloid is located along the y-axis 
but displaced by j/o from the origin as shown in Fig. ^ 
The surface of the paraboloid can be written as: 



(2) 

We define the bow shock region to 



y = QyP + yo 

where P — ^ 
be the volume between this surface and the same surface 
displaced by d (thickness) in -|-y direction. 

The flow of the gas in the CW zone is assumed to be 
tangential to the surface of the paraboloid; hence, there is 
no azimuthal component. Then, the velocity (vcw) of the 
gas flow at point P in the CW zone can be written as: 

Vcw=Wj,y + u;l (3) 

where 1 = cos (pci + sin 0cX. Vy and vi are the velocity 
components which are parallel and perpendicular to the y 



^ See [Huang fc Weigert (1982|) , [Cirard fc Willson (1987|) and 
[Canto etaL(1996 ) for the predicted geometry of the thin-shell 



shocks created by the interaction of the two stellar winds in a 
highly radiative binary. 
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Fig. 1. This figure illustrates the model configurations 
of V444 Cyg. The W-R star is placed at the center of 
the cubic boundary, and it is surrounded by a spherical 
atmosphere which is consistent with the CMFGEN spec- 
tral model. The O star is located below the W-R star 
in the diagram, and is covered by the tilted paraboloid 
bow shock with a given thickness. The strong stellar wind 
from the W-R star is interrupted by the shock front, and 
the density behind the shock is assumed to be negligibly 
small for simplicity. The directions of an observer at phase 
= 0, 0.25, 0.5, 0.75 are indicated at the top, bottom, left 
and right edges since the orbital incUnation is about 80° 
(almost edge-on view). 



axis, respectively, ay is a parameter to control the asymp- 
totic open angle of the bow shock. Further, by taking the 
derivative of Eq. ^ and using the absence of an azimuthal 
component of Vcw, we obtain: 



dl vi 



(4) 



If the magnitude of the velocity (vcw) at point P is 
similar to that of the spherical flow (vr) of the pre-shock 
W-R wind at P , i.e., Wcw ~ Vr, then 



2 2 I 2 



(5) 



Since the flow of gas in the bow shock region is usually 
slower than the surrounding flow according to the hydro- 
dynamics model of Pittard fc Stevens (1999|) , we introduce 



a free parameter, and rewrite the equation above as: 



S Vr 



,1/2 




Fig. 2. This illustrates the simple geometric model of the 
bow shock and the flow of the gas. The W-R star is located 
at the origin (0), and the O star is located just behind 
the bow shock on the y axis. The gas in the bow shock 
region is assumed to flow tangentially to the surface of the 
paraboloid in the +y direction. The gas velocity v^w at a 
point p on the surface is decomposed into Vy and u/, but 
there is no azimuthal component (i.e., no rotation about 
the y axis). 



In general s < 1, and is a function of the position (e.g., 
s — s (y)). However, in most of the analysis, s — 1/2 is 
assumed for simplicity. 

Using Eqs. ^ ^ and the components of Vcw can be 
expressed, in terms of the position (y, I) on the paraboloid, 
as 



2ay I s Vr 



{2ayiy 



1/2 ' 



and 



vi 



l + i2ayiy 



1/2 ■ 



(7) 



(8) 



These expressions are then corrected for the tilt of the 
bow shock due to the fast orbital motion of the system, 
as mentioned earlier. A more accurate velocity field in the 
bow shock region is needed to properly model the line 
profile and the line polarization, but it is not important 
for the continuum calculations. 

Considering the radial dependence of the density, we 
define: 



(9) 
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where r is the distance between the origin (0) and a point 
p on the surface of the parabola, and yo is the distance be- 
tween the origin and the head of the paraboloid as shown 
in Fig. ^. The index n is a scaling parameter, and po is the 
density at the head of the parabola (at y = yo). This for- 
mulation is similar to one in [Stevens fc Howarth (l"999| ). 



2.5 



Po is estimated from the amount of the excess emission 
seen in the He i 5876 A line, and is discussed in the later 
section (§ 2.2.3). In our analysis, we use the bow shock 
density expressed in Eq. ^ with n = 2, and the thickness 



(d) of the bow s hock w 0.1 i?o, estimated from Pittard 



& Ste\ ens (1999 ). In predicting the contimmm polariza- 
tion level, the global distribution of the electron density is 
more important than small structure. The column density 
or the electron scattering opacity through the bow shock 
is related to the amount of gas in the cavity which will be 
absent if there is no bow shock. 

Some sensitivity of the line strengths and polarization 
may arise from assuming whether the helium is singly 
or doubly ionized. Hei 5876 A and Hen 5412 A lines 
are often referred to as "diagnostic lines" for model- 
ing a WN star atmosphere (e.g., Hamann et al. 1995| ). 
The latter is an optically thick line; hence, the absorp- 
tion and scattering terms in radiative transfer equations 
must be treated properly. The gas in the bow shock is 
ra,pid]y cooled by radiative processes, a.nd it could 



tentially cause a significa.nt amount of absorption in the 
He I 5876 A line. Figure. ^ shows the approximate loca- 
tion where Hei 5876 A and Heii 5412 A lines originate ac- 
cording to a CMFGEN model of a typical WN5 star with 
RwR = 5 i?Q . The figure shows that the He ii emission 
peaks at around 3 i?wR from the center of the W-R star, 
and the He i emission peaks around the location where the 
O star resides. This explains why the emission from the 
CW zone is more clearly seen in the He i line. In the fol- 
lowing, the variability of the He i line due to the presence 
of the bow shock will be examined. 



2.1.1. Effect of the Colliding Winds 

As explained earlier, the emission from the bow shock re- 
gion is prominent for some lines (s 



Marchenko et al. 



1994). In this section, the expected effect of the bow shock 



emission on the He i 5876 A line will be demonstrated us- 
ing the 3-D Monte Carlo model. He i 5876 A is primarily 
a recombination line; hence, to first order we assume the 
line emissivity in the bow shock region is simply propor- 
tional to the square of the density. The tilt angle S — 15° 
and the open angle = 90° (Marchenko et al. 1994) are 
used for the geometry of the bow shock. The left graph 
in Fig. ^ shows a sequence of He I 5876 A emission lines 
as a function of phase (0 — 1 with steps of 0.1 from the 
bottom) arising from a V444 Cyg model which includes 
an unrealistically strong bow shock emission. In fact, the 
line emission is dominated by the bow shock emission, and 
the underlying line emission from the W-R atmosphere is 
barely visible in this figure. This plot demonstrates the 
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Fig. 3. Illustration of the approximate location where 
Hei 5876 A (solid) and Heii 5412 A (dashed) lines orig- 
inate according to a CMFGEN model of a typical WN5 
star with i?wR = 5 Rq- The amount of line emission orig- 
inating in the interval dlog(r/i?wR) is proportional to 
^ dlog (r/i?wR) where r is the distance from the center of 
the W-R star. The solid vertical line at log(r/i?wR) = 
0.88 indicates the location of the O star assuming that 



the separation of system is a = 38 Rq ( Marchenko ct al 



1994). According to this figure, the Heii emission peaks 



at around 3 i?wR from the center of the W-R star, and 
the He i emission peaks around the location where the O 
star resides. 



basic behavior of the emission from the bow shock region. 



(See also Liihrs 1997 and Bartzakos et al. 2001 for similar 
models and discussion.) Ai (f) — 0.8 — 1.0, the bow shock 
is pointing away from an observer and tilted to the light 
of sight (Fig. 1^); therefore, the emission peak appears on 
the red side of the flat-top He i emission line. The shape 
of the CW emission changes very rapidly from </> = to 
0.1 because a rapid change in the velocity component of 
the bow shock arms along the line of sight occurs near 
(j) ~ 0.1. The emission from the bow shock flattens out 
during = 0.1 — 0.4 and (j) = 0.7 — 0.8 since the bow shock 
arms are almost perpendicular to the direction of the ob- 
server. At = 0.4 — 0.5, the peak appears on the blue 
side since the bow shock arms are pointing towards an ob- 
server. A sequence of He i 5876 A line profiles with a more 
realistic level of bow shock emission is shown in the right 
panel of Fig. ^. The figure shows that the line strength 
increases around 4) — and 0.5 because the continuum 
level drops during the primary and secondary eclipses. 

2.1.2. star Absorption Effect 

The observed spectrum of V444 Cyg shows absorption 
lines due to the O star atmosphere on top of the broad 
emission lines from the W-R atmosphere. The position of 
absorption shifts as the relative speed of the stars changes. 
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Fig. 4. Left: The phase dependent emission ft'om the bow 
shock region on top of Hei 5876 A hne. The emissivity 
of the bow shock is assigned to be unreahstically strong 
in order to demonstrate its effect more clearly. Right: As 
for the left figure, but the emissivity of the bow shock is 
reduced by a factor of ~ 10. This figure qualitatively de- 
scribes the variability seen in the observation of this line 

The profiles are plot- 
1 (top) in 0.1 



(e.g., see Marchenko et al. 1994). 
ted from phase = (bottom) to phase 
phase steps. The profiles are normalized to continuum, and 
shifted upward by 2 (right) and 0.12 (left) consecutively 
at each phase. In these models, the tilt angle (5 = 15° and 
open angle = 90° ( Marchenko et al. 1994[ ) are used. 



The distortion caused by this effect is prominent for some 
lines at certain binary phases. Since the model does not 
take into account the photospheric spectrum of the O 
star, the O star absorption lines have been added to the 
CMFGEN model in an approximate manner. In principal 
they could easily be incorporated into the Monte Carlo 
method by changing the appropriate boundary conditions. 

The O star component of V444 Cyg has spectral type 



05-06.5 III-V (Marchenko et al. 1994). The O star (06) 



spectrum is modeled separately by CMFGEN, and the 
effect of the rotational broadening of lines is added to the 
spectrum. In this process, the rotational s peed Vrot sin i = 
215(±13)kms-i ([Marchenko et al. 1994| ) is used. 

Fig. H shows the combined spectrum of the W-R and 
the O star computed by CMFGEN as a function of the 
orbital phase. The sinusoidal motion of the O star absorp- 
tion component is clearly seen in both lines. To examine 
only the effect of the O star, the emission from the CW 
zone is not included in the model spectra shown in this 
figure. A specific luminosity ratio {q = L-yvR/io = 0.45 
at A = 5630 A) was used when the two spectra were com- 
bined. 




5400 5450 
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5850 5900 
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5950 



Fig. 5. The figure illustrates how the O star absorp- 
tion affects the appearance of the Hei 5876 A and 
He II 5412 A lines from the W-R atmosphere as a function 
of the orbital phase {(f)) . The spectrum of the W-R and the 
O stars computed with the spherical model, CMFGEN, 
are combined using the monochromatic luminosity ratio, 
q = Lwr/^o = 0.45 at A = 5630 A. The figure does not 
include the effect of the colliding wind for clarity. The pro- 
files are plotted from (j) = (bottom) to <j) = 1 (top) with 
0.1 phase steps. 

2.1.3. Combined Effects Compared with Observation 

In order to use CMFGEN for the model fitting, we must 
choose a binary phase when the lines are least affected by 
the bow shock emission and the O star absorption. Fig. || 
shows that the absorption by the O star does not affect the 
overall shape of Hei 5876 A, but it does affect the shape 
of He II 5412 A. At w 0.2 and (f) « 0.7, the Hei line 
shape is least affected by the O star, but at these phases 
the emission from the bow shock region is spread over the 
Hei line as seen in Fig. ^. As a result, the line flux level 
is likely slightly higher than that of the spherical model. 
On the other hand, at ~ and (j) « 0.5, the emission 
from the bow shock region is well concentrated near the 
blue or the red edge of the He i line, but at these phases, 
the stars are eclipsed; hence, the line strength will be 
underestimated in a simple spherical model (CMFGEN). 
Meanwhile at « 0.3 and (j) « 0.8, the bow shock emis- 
sion is concentrated only on either the blue or the red 
side of the Hei line, and the O star absorption affects 
the same half side of the line profile. Unfortunately, at 
(j) « 0.3, the O star absorption is on the red side of the 
helium line, and the bow shock emission is on the blue 
side of the line. So the line is contaminated on both sides. 
On the other hand at (/> « 0.8, both the O star absorption 
and the emission from CW are on the red side of a line; 
hence we can reliably use the blue part of He i 5876 A and 
Hell 5412 A for the model fitting. Figure |6| shows the 
sample CMFGEN model spectrum compared with the ob- 
served spectrum of Marchenko (priv. comm) at six differ- 



8 



Kurosawa et al.: Polarization and Mass-Loss Rate of V444 Cyg 



ent phases. The agreement between the model and the ob- 
servation at w 0.8 is excellent. Near the eclipsing phases 
(0 « and 0.5), the spherical model CMFGEN does not 
fit the observation well since the contribution from the 
bow shock is not included in the model. The two narrow 
Na I interstellar absorption lines are present in the red side 
of observed Hei 5876 A lines. 



2.2. Mass-Loss Rate Estimation of the W-R 
Component 

The idea behind the mass-loss estimation performed here 
is very similar to the one in STLl and STL2. We assume 
that there is a correspondence between the amount of con- 
tinuum polarization and the mass-loss rate of the W-R 
component. This is supported by the analytic expression 
of M in terms of Ap (a semi-major axis of polarization on 
Q-U plane) derived by STLl, using the basi c results from 
the classic paper on binary polarization by Brown et al 
l(197l ). 



2.2.1. Investigation of the W-R Radius 

For an assumed mass-loss rate, M, of the W-R star, there 
are several CMFGEN models which can fit the observa- 
tional spectrum of Heii 5412 A and Hei 5876 A simul- 
taneously. As discussed in the previous section, the spec- 
trum at (/) « 0.8 is used in the model fitting since the 
contamination from the O star atmospheric absorption 
and the CW emission is expected to be small at this or- 
bital phase. Fig. shows the fits of the observed spec- 
trum at = 0.826 for different values of i?wR (=2-5, 4.0, 
5.0 i?Q)Q with the mass-loss rate of the W-R star fixed 
(MwR = 0.6 X 10"^ Mq yr^^). The corresponding volume 
filling factors of the models are 0.08, 0.04, 0.05. All three 
models fit the observation very well; therefore, the spec- 
tral solutions are degenerate. The He ii profiles with differ- 
ent values of RwB. a-re not significantly different from each 
other when the O star spectrum is added to the W-R spec- 
trum with an appropriate value of q (= 0.30, 0.35, 0.45 
for Rys ii — 2.5, 4.0, 5.0/^0 models respectively). In the 
case or a single W-R star spectrum, a small but noticc- 
able difference in the profile shape of the Hen line, usu- 
ally becoming more triangular for a larger i?wRj should 
be present. Unfortunately, it is hard to see this effect in 
the V444 Cyg model because the continuum flux of the 
O star component is stronger {q « 0.5); hence, the line 
strength is weakened. 

In an attempt to remove this degeneracy, the light 
curves (/, Q, U) are computed by the 3-D Monte-Carlo 
model (Kurosawa & Hillier 2001) for the spectral models 
shown in Fig. |7| for fixed i?o and Mwr values. First, mod- 
els without the bow shock (Figs. || - ^2|) are considered 
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The different solutions correspond to different assumed dis- 
tances for V444 Cyg. The W-R radius (i?wR) used in this paper 
corresponds to the inner boundary of the model atmosphere 
where the Rosseland optical depth tr ~ 30. 



Fig. 7. Left: the model fits of Hen 5412 A for three dif- 
ferent radu (2.5, 4.0 and 5.0 i?©) of the W-R star. Right: 
same as the left, but for Hei 5876 A. The profile shapes 
of both lines fit the observation well for a relatively wide 
range of the W-R star radius. When the O star spectrum is 
added to the W-R spectrum with an appropriate value of 
q, the profiles with different values of i?wR are not signif- 
icantly different from each other. The values of q used are 
0.30, 0.35 and 0.45 for i?wR =2.5, 4.0 and 5.0% models 
respectively. The corresponding volume filling factors of 
the models are 0.08, 0.04, 0.05. 



for clarity. The effect of the bow shock on the light curve 
will be examined later. A small difference is expected to 
be seen (c.f., STL2) in the shape of the polarization curve 
near the secondary eclipse (0 « 0.5, i.e., when the W-R 
star is behind the O star). 

The results, using Rq = 7.2 Rr^ and A/wr = 0.6 x 
10~^MQyr~^, are shown in Fig. ^. The orbital inclina- 
tion angles used for the fits are i = 78.8°, 78.5° and 
78.0° for the models with i?wR = 2.5 i?0, 4.0 i?0 and 
5.5 i?0 respectively. The inclination angles are consistent 
with i = 78° ± 1° dChercpashchuk 1975| ) obtained from 
the light curve solution. STL2 estimated the inclination 
angle by fitting the double sine curves of Q and U po- 
larization components with the formula given by Brown 
et al. (1978"! ). Their resuh is i = 80.8° ±1.6°. In Fig. 
the three models fit the observed / light curve moder- 
ately well. Although not shown here, all three models also 
fit well the observed / light curve outside of the phase 
range shown in Fig. ^. The deviations from perfect anti- 
symmetry around (j) — Q.b seen in the observational / light 
curve and the anomaly seen around (j) « 0.53 — 0.54 in the 
U observational light curve are not well fitted, and may 
be related to the presence of the bow shock in the system. 

The width of the / light curve around (j) = 0.5 is 
slightly wider for the model with the larger W-R star ra- 
dius. The i?wR = 5.0 i?Q model fits the overall shape 
of the light curve better than the two other models with 
smaller radii. However, the width of the secondary eclipse 
depends also on the mass-loss rate of the W-R star, and 
we discuss this dependency later. Unfortunately, no sig- 
nificant difference in the Q and U light curves are seen for 
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Fig. 6. The figure shows one of the best fit CMFGEN models (Dotted) at 6 different binary phases. The observational 
data (Solid) were provided by Marchenko (priv. comm.). The O star model spectrum (06V) is added to include the 
effect of the atmospheric absorption effect. The monochromatic luminosity ratio q = L^^^/Lq = 0.45 at A = 5630 
A is used when the W-R star model spectrum is combined with the O star model spectrum. In general, the fits are 
good at all phases shown, especially at « 0.3 and (f> « 0.8. Two narrow Nai interstellar absorption lines are present 
in the red side of observed Hei 5876 A lines. The emission from the bow shock is not included in the model shown 
here. 



different values of W-R star radii, unlike that found by 
STL2. In fact, the difference in Q and U light curves seen 
by STL2 for different i?wR (see their Fig. 8) is not purely 
due to changing Rwr, being also coupled with Af-yvR and 
the adopted inclination. 

In this analysis, we were not able to constrain the value 
of Rwr- Models with a wide range of Rwr = 2.5 — 5.0 Rq 
can fit the observed /, Q and U light curves reasonably 
well, and there is no significant difference among them as 
long as appropriate values of the orbital inclination angle 
and the O star radius are chosen for each model. Despite 
this uncertainty in the W-R star radius, we continue to 
derive the mass-loss rate of the W-R star. 

2.2.2. Determination of the W-R Mass-Loss Rate 

Next, the observed helium spectrum is modeled for dif- 
ferent values of Mwr by separately adjusting the values 



of / and q for Rwr = 5.0 and 2.5 i?Q. Fig. || shows the 
results for Mwr = 0.3, 0.6 and 1.2 x lO'^MQyr^^ with 
the fixed radius of the W-R star: i?wR = 5.0 i?Q. No sig- 
nificant difference in the profile shapes is seen in these 
models. To remove this degeneracy, once again, the po- 
larization curves were computed; however, this time the 
amplitude of the "double-sine" shaped polarization curves 
were fitted. The polarization curves for each value of Mwr 
corresponding to the spectral models in Fig. |^ are shown 
in Fig. |l^ along with the light curves at A = 5630 A. 
The model parameters used for the fits in Fig. |l^ are 
summarized in Table ||. The figure clearly shows that the 
model with Mwr. — 1.2 x 10~^MQyr~^ overestimates 
and the one with AfwR = 0.3 x 10~^MQyr~^ underes- 
timates the amplitude of the observed Q and U curves, 
while Mwr = 0.6 x 10^^ Mq yr^^ model is consistent 
with the observations. The width of the / light curve for 
A — 5630 A around the primary eclipse (</) = 0.0) with 
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1.05 



? 0.95 - 




Fig. 8. The continuum flux, Q and U polarization near 
the secondary eclipse (where the O star is in front of the 
W -R star) are computed by the 3-D Monte-Carlo model 
of Kurosawa fc Hillier (200l| ) for the three spectral models 
shown in Fig. |7|. Circles: Observation, Dot: i?wR — 2.5 i?Q 
hq = 0.30 model, Dash-Dot: i?wR 4.Oi?0 & <? = 0.35 
model, Sohd: i?wR = 5.0 i?0 & g = 0.45 model. The typ- 
ical size of the error in model calculations of Q and \J is 
indicated near the upper left corner of the middle and the 
bottom plots. The optical light curve data and optical po- 
larization data (V band) are from Kron & Gordon (1943) 
and STL2 respectively. 



-^^WR = 0.3 X 10 ^ Mq yr ^ is slightly narrower than the 
observation, and that with MwR = 1-2 x 10~^MQyr^^ 
is slightly wider than the observation. The observed / 
light curve at A = 5630 A is fitted with the Mwr = 
0.6 X 10~^ Mq yr~^ model very well. 

Since there is a large uncertainty in the -RwR value, 
the results of similar calculations with i?wR = 2.5 i?Q are 
shown in Figs. |l^ and 12, They are very similar to the re- 
sults of the i?wR = 5.0 i?Q models, and again the model 
with A%R = 0.6 X lO-^M0yr-i fits the observed /, Q, 
and \J light curves at the wavelength around 5630 A con- 
sistently. A summary of the model parameters is given in 
Table ll 



2.2.3. Effects of the Bow Shock on Light Curves 

We now consider the possible eff'ects of the presence of the 
bow shock region on our computed light curves (Figs, 
and p^ ). To do so, the observed Hei 5876 A line at the 
phase where the effect of the CW emission is most promi- 
nent is modeled first. The spectrum taken at </) = 0.419 
(see Fi g. ^ is chosen for this purpose. The profile at this 
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Fig. 9. The observed Heii 5412 A (left) and 
He I 5876 A (right) are fitted with the models for different 
values of the W-R star mass-loss rate. i?wR = 5.0 i?Q is 
used for all the models shown here. Solid: observation. 
Dot: Model A (Mwr = 0.3 x lO-^MQyr"!), Dash: 
Model B (Mwr = 0.6 x lO'^MQyr^^), Long Dash: 
Model C (Mwr = 1-2 x lO-^MQyr"!). See Table. | 
for other parameters used. The profiles from all three 
models are very similar to each other, and their fits to the 
observation are reasonable. It is impossible to distinguish 
between the three models. 
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phase js fitted by adjusting po in Eg. ^ to produce a sim- 
ilar amount of the excess emission as seen on the top of 



Fig. 11. Same as Fig. but with i?wR — 2.5Rq. Solid: 

observation. Dot: Model D (Mwr = 0.3 x 10"^ Mq yr'^), 
Dash: Model E (Mwr = 0.6 x 10"^ Mq yr^^), Long-Dash: 
Model F (A/wR = 1.2 X lO^^MQyr-i). See Table | for 
other parameters used. 



the He I 5876 A line. Since the true line opacity and emis- 
sivity in the bow shock region are not known, this was 
done just to produce the bow shock line emission seen in 
He I 5876 A. We do not intend to perform a very rigorous 
fit of this line although that would be a good topic for a 
future investigation. The result from the 3-D Monte Carlo 
calculation for Model B is shown in Fig. |l^. A tilt of the 
bow shock, 6 — 15°, and open angle = 90° ( Marchenko 
et al. 1994| ) were used in this model. The fit shown in 



Fig. |l3|is reasonably good except for the flux level around 
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Table 2. Model Summary 
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MODEL 


A 


B 


C 


D 


E 


F 


RwR 




5.0 


5.0 


5.0 


2.5 


2.5 


2.5 


LwR 




2.0 


2.0 


2.0 


2.0 


2.0 


2.0 


MwR 


[xlO-^Meyr-i] 


0.3 


0.6 


1.2 


0.3 


0.6 


1.2 




1700 


1700 


1700 


1700 


1700 


1700 


/ 




0.018 


0.05 


0.15 


0.025 


0.075 


0.225 


q (A5630) 


0.52 


0.45 


0.40 


0.31 


0.31 


0.29 


q (A22000) 


1.15 


1.10 


1.05 


1.13 


0.96 


0.90 


Ro [Rq] 


5.7 


7.2 


10.0 


5.5 


6.9 


8.6 


i [deg 


] 


80.5 


78.0 


73.0 


82.0 


79.5 
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Fig. 10. For each model shown in Fig. ||, the polarization and flux variations are computed. The three plots on the left 
are the relative flux, Q and U polarization at A = 5630A as a function of the binary phase. The three plots on the right 
side are the same as those on the left, but computed at A = 2.2 /im. Circle: observation. Dash Dot: Model A (Mwr = 
0.3 X lO^^MQyr-i), Solid: Model B (AfwR = 0.6 x 10"^ A/q yr^^). Dot: Model C (AfwR = 1.2 x lO^^^ yr'^). 
The amplitude of the observed polarization curves is not consistent with either Model A or Model C. On the other 
hand, Model B (AfwR = 0. 6 x lO'^Af^yr"^) fits the observation very well. The optical light curve data and optical 
polarization data are from Kron & Gordon (1943) and STL2 respectively. The observed infrared light curve is from 
Hartmann (1978|). 
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Fig. 12. Same as Fig. |T^, but the models are calculated for Rwr — 2.5 Rq. The spectral models corresponding to 
these models are shown in Fig. |ll]. Circle: observation, Dash-Dot: Model D (MwR = 0.3 x 10~^ Mq yr~^), Solid: Model 
E (MwR = 0.6 X 10"^ Mq yr-i), Dot: Model F (AfwR = 1-2 x lO"'^ yj.-i). Model D (Mwr = 1.2 x lO"'^ yj.-^) 
and Model F (Mwr = 0.3 x 10~^MQyr~^) do not fit the amplitudes of the observed polarization curves. On the 
other hand, Model E (Mwr = 0.6 x 10~^ Mq yr~^) fits the observation very well. 



V — 200 — 500 km s^^ where the observed spectrum is af- 
fected by the O star absorption. For simplicity, O star 
atmospheric absorption effects are not included in this 
model; however, an appropriate q was used to produce the 
correct O star continuum flux. The ratio the bow shock 
density to the background density (pcw/pbg), due to the 
W-R stellar wind, is found to be about 40 according to this 
model. This ratio is similar to that of Pittard & Stevens 



(19990 



With this bow shock model, the light curves (/, Q, 
and U) at X = 5630 A of Model B (Mwr = 0.6 x 
10~^ Mq yr~^ and i?wR — 5.0 Rq) in Fig. |l^ are recom- 
puted. Fig. |lj shows the result. There is no significant 
effect of the bow shock region on the continuum /, Q and 
U light curves. The two models are identical to each other 
within the range of the error. The same conclusion was 
obtained for the other models shown in Figs. and n2. 



2.2.4. Effect of the 0-Star Wind using a 
Hydrodynamical Model 

To check the validity of ignoring the O star wind and 
the geometrical bow shock model, the polarization light 
cu rves were computed using the hydrodynamical model 
of Pittard fc Stevens (1999 ), which includes the orbital 
motion of the stars and a parameterized form of ra- 
diative driving. First, the opacity of the W-R star at- 
mosphere was comp uted with CMFGEN usin g the pa- 
rameters adapted by Pittard fc Stevens (1999 ): -Rwr = 



4.0 i?0, Ro = 10.0 i?0, L 



^WR 



1.40 X IO^Xq, Lc 



2.4 X 105^0, MwR = 9.72 x IQ-^ M^yr-^ Mq = 
5.78 x 10-'^ Mqyt-\ T/WR = 1800kms~i and = 
2900kms-i. The CMFGEN model with these parameters 
does not fit the observed helium spectrum, but the opacity 
from this model will be used in the polarization calcula- 
tion to maintain the consistency with the hydrodynamical 
model. 

Using this W-R star atmosphere model, the following 
two models are computed: l.a model without the bow 
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Fig. 13. Spectral model of Hei 5876 A at phase = 0.419 
where the emission from the bow shock is prominent. The 
profile is fitted by adjusting po in Eq. ^to produce a simi- 
lar amount of the excess emission as seen on the top of the 
He I 5876 A line. T he tih of the bow shoc k 5 = 15° and 
open angle = 90° ( Marchenko et al. 1994 ) were used in 
the geometrical bow shock model introduced in § 2A_. To 



compute the profile, the emissivity and the opacity from 
Model B were used in the 3-D Monte Carlo model. 
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Fig. 14. Comparisons of the light curves of Model B in 
Table ^ with and without the bow shock. The opac- 
ity of the bow shock used here is obtained by fitting 
He I 5876 A at phase = 0.419 as shown in Fig. |l|. Solid 
line: Model B without the bow shock. Dot-Dash: Model 
with the bow shock. The data points are indicated by cir- 
cles. Within the range of the error, the two light curve 
models are identical to each other. The size of the typi- 
cal error (arising from random fluctuations) for the model 
Q and U polarization curves are shown on the upper left 
hand corner of each plot. The size of the error for / is too 
small to be shown on the plot. 



shock and the O star wind, 2. a model with the O star 
wind and the bow shock from the hydrodynamical calcu- 
lation of Pittard & Stevens (1999). The density distribu- 
tion for the second case is shown in Fig. |l^. To assign the 
opacity and emissivity in the O star wind and the bow 
shock regions, we assumed that the thermal opacity and 
the emissivity are proportional to the square of the den- 
sity, and the electron scattering opacity is proportional 
to the density. Then the opacity and the emissivity are 
scaled with those from the W-R star atmosphere model. 
The light curves computed for these models are shown in 
Fig. |l6|. No significant differences between the two mod- 
els are seen in this figure. This validates one of our main 
model assumptions, i.e., the O star wind does not signif- 
icantly affect the continuum polarization curves. It also 
assures that the details of the bow shock structure are 
not important for predicting the continuum polarization 
curves though they are very important for a line calcula- 
tion. 

The model light curves in Fig. |l^ do not fit the ob- 
served light curves. By increasing the inclination angle, 
the depth of the eclipses (at ^ = and 0.5) will increase; 
however, without adjusting other parameters, it is impos- 
sible to decrease the depth of the secondary eclipse. It 
is clearly seen from the location of minima in the eclips- 
ing Q and U light curves {(j) ~ 0.46) that the O star ra- 
dius is too large. Surprisingly, the amplitudes of the Q 
and U double sine curves fit the observation well with 
Mw,j^ = 9.72 X lO"'^M0yr"^ This seems contradictory 
to our earlier analysis, but is due to the inconsistent W-R 
model, which does not fit the observed helium spectrum, 
being used in the polarization calculations. Since the pa- 



rameters from the model of Pittard & Stevens (1999) are 
fixed, q — Lwr/Lq — 0.583 is also fixed in these models. 
li q = 0.38 is used by increasing Lq instead, the spec- 
tral model will fit the observed helium spectrum, but the 
higher Lq results in larger amplitudes of the double-sine 
polarization curves which then no longer fit the observa- 
tion. 

In summary, from this analysis, we conclude that the 
mass-loss rate of the W-R star in V444 Cyg is 0.6 (±0.2) x 
lO^^MQyr"^. The size of the error in M is estimated 
from the light curves in Figs. |l^ and|l^. The mass- loss rate 
determined here is consistent with the values of STL2, and 
it is slightly higher than that of Underhill et al. (19901 ). 



(See Table |^ for other published values.) The range of the 
volume filling factor for the W-R star atmosphere is esti- 
mated to be 0.050 — 0.075 with the corresponding range 
of the W-R star radius, 5.0 — 2.5 i?Q. We also found that 
the presence of the bow shock, the exact details of the 
bow shock, and the presence of the O star wind are unim- 
portant for the continuum fits, and they do not affect the 
mass-loss rate estimate given here. 



3. Discussion 
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Fig. 15. This figure shows the density distribution of the 
V444 Cyg model described in § 2.2.4. The density, includ- 



ing the bow shock, near the O star (upper-left) is from the 
hydrodynamical model of Pittard fc Stevens (1999| ), and 



is added to the spherical density distribution of the W-R 
star (center) derived from a CMFGEN model. The param- 
eters used in this model are: i?wR = 4i?Q, Rq = 10 Rq 
anda = S8Rq. Lwr = 1.40 x 10^ Lq, = 2.4x10^2.0, 
MwR = 9.72xlO-^MQyr-\Mo = 5.78x10~^Mq yr'^, 
p-WR ^ i800kms-i and V2 = 2900kms-i. 



Fig. 16. The light curves are computed with the parame- 
ters used in the hydrodynamical calculation of Pittard & 
Stevens (1999| ). They are: i?wR = 4.0 i?0, Rq = 10.0%, 
Lwr = 1.40 x IO^Xq, Lq = 2.4 x IO^Lq, Mwr = 
9.72 X lO-^Mgyr-i and Mq = 5.78 x lO-'^Mgyr"!. 
The solid lines are from the model without the bow shock 
and the O star wind. The dashed lines are the model with 
the O star wind and the bow shock from the hydrody- 
namical model of [Pittard fc Stevens (1999| ) (Fig. |l|). No 
significant differences between the two models are seen in 
this figure. 



the broadband (K) photometry obtained by Hartmann 



3.1. I mrared Light Curves 



(1978). The depth of the primary eclipse (at = or 



In order to check the consistency of our best fit models 
(Models B and E with Mwr = 0.6 x lO'^ MQyr~^) in 
§ |2.2| with some other observational aspects, the infrared 
light curves (at A — 2.2 ^m) were also computed. In or- 
der to calculate the IR light curves, the monochromatic 
luminosity ratio q ~ Lwr/Lq at A — 2.2 /im must first 



be esti nated. In Model B, q was determined to be 0.45 at 
A — 5630 A in the previous section. By normalizing the 
O star fiux with this value of g at A = 5630 A, q for the 



1 when the W-R star is in front of the O star) is fitted 
well for both Rwr = 2.5 Rq and 5.0 i?Q models, but 
the secondary eclipse modeled (at </) = 0.5 when the O 
star is in front of the W-R star) is shallower than the 
observed IR light curve. We also notice that the maxi- 
mum I flux level (/max) in the observed light curve is not 
well defined. Cherepashchuk ct al. (1984) and Hamann 
& Schwarz (1992) had a similar problem - the primary 



and secondary eclipses are too shallow in their light curve 



wav elefigths between 3000 A and 33, 000 A is plotted in 
Fig. |l7| for three different O star effective temperatures. 
The figure shows that the O star is brighter than the W- 
R star at A = 5630 A, but the W-R star is brighter than 
the O star after A ~ 20,000 A. In model B, an O star 
temperature of Tq — 35, 000 K was used; hence, the cor- 
responding value of g at A = 2.2 fim is 1.10 according to 
the graph. For models with higher O star temperatures, 
the values of q are slightly higher. 

The resulting model IR light curves are shown in the 
right hand side of Figs. 10 and |lj along with Q and U po- 
larization light curves for M-yvR = 0.3 x 10~^ Mq, Mwr = 
0.6 X 10-^ Mqjt-^ and Mwr = 1.2 x lO"'"^ Mq yr^^. 
The figure also includes the observed IR light curve from 



model. While the light curve model of Hamann & Schwarz 
(1992D did not fit the depth of both eclipses, our model is 



consistent with at least the primary eclipse. 

Cherepashchuk et al. (1984) discussed two possible 
causes for this problem: a), the existence of "third light" 
and b). the result of ignoring the reflection effect. They 
noted that one possible source of third light is the contri- 
bution of emission lines in the broad-band photometry at 
2.2 /xm as originally mentioned by Hartmann (1978| ), but 
we now expect this effect to be generally small (unless the 
He I 2.06 iim emission line is abnormally large). To exam- 
ine this possibility, we have run a model which includes 
the reflection by the disk of the O star. However we find 
that it does not affect the in frared and optical light curves. 
[Cherepashchuk et al. (1984 ) also conjectured that the in- 
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Fig. 17. The continumn flux ratio {q = Fwr/Fq) of 
Model B as a function of wavelength for three different 
O star I effective temperatures. The lines were normalized 
at A = 5630A where the luminosity ratio q — 0.45 was 
determined from the spectroscopic model fit (Fig. 



crease in the depth and the width of the secondary eclipse 
(from optical to IR) is a result of clumps in the WR stellar 
wind. 

If the wrong flux ratio at A = 2.2 /im is used in our 
model, the secondary eclipse could become shallower at 
this wavelength. The differences in the depth of the sec- 
ondary eclipse in Figs. |l^ and 12 are approximately 8% 
and 5%. This could be explained by the wrong O star 
temperature used in the model. According to Fig. |l^, q at 
2.2 yLtm also increases about 8% if we use Tq = 45, 000 K 
although this would not fit the absorption components of 
the observed spectrum. 





Model B 


Model E 


7?w« li?^l 


5.0 


2.5 


no [Hq\ 


7 


f< Q 


AiwR IXlU AIq yr 1 


0.6 


0.6 


/ 


0.05 


0.075 


q (A5630) 


0.45 


0.31 


q (A22000) 


1.10 


0.96 


i [deg.] 


78.0 


79.5 


(O star) 


-5.20 


-4.85 


My (W-R star) 


-4.33 


-3.58 


My (O-HW-R) 


-5.60 


-5.15 


m„ (0-hW-R)(*) 


8.38 


8.84 


DM (**) 


11.1 


10.6 





0.68 Lundstrom 



(*) Assuming a dist ance of 1.7 kpc and Et^ 
& Stenholm (1984) - They assumed the intrinsic color to be 
(6 -»;)„ = -0.30. 

(**) Using the observed value of rr iy — 8.27 and Ay — 2.79 
from Lundstrom & Stenholm (1984). 



3.2. Polarization Curves near the Secondary Eclipse 

In both Model B (Fig. [l^) and Model E (Fig. |l|), the am- 
plitude of the U polarization variation near the secondary 
eclipse {</) = 0.5) underestimates the observed level. A 
larger polarization amplitude can be obtained if the ra- 
dius of the O star is increased, although this will result in 
an inconsistent depth of the secondary eclipse. A possible 
cause of this inconsistency is the inaccurate stellar wind 
structure very close to the core of the W-R star. For all 
the models in Ta ble p|, /3i = 1 was use d in the modified 
beta- velocity law ( Hilher fc Miller 1999| ) (see Eq. ||) where 
Pi is equivalent to P in the classic beta-velocity law for 
the "inner" part of the stellar wind. To see if changing the 
structure of the inner W-R stellar wind affects the com- 
puted continuum polarization level, a model with /3i = 2 
(a slower wind) is also calculated. Unfortunately, the am- 
plitude of the polarization curves (not shown here) near 
(f) = 0.5 did not change significantly for the new Pi value. 



Changing the value of Pi also did not influence the depth 
of the secondary eclipse in the IR light curve. 

We note that STL2 were able to fit the observed Q and 
U variations around (f) = 0.5. There are two reasons for 
their success: 1. they did not fit the observed / light curve 
consistently, and 2. they use a fixed orbital inclination an- 
gle. In their fitting procedure, the radius of the O star 
was used as a fitting parameter, but changing the O star 
radius changes the depth of the eclipse around 4> = 0.5. 
Consequently, the inclination angle must also be adjusted 
when the O star radius is changed, to maintain a consis- 
tency with the observed / light curve. We should also note 
that the amplitudes of polarization curves around 4> = 0.5 
are very sensitive to the orbital inclination. STL2 adopted 
i = 78.9° from Robert et al. (199C ). However, if instead 
i — 78° from Cherepashchuk (1975) were used, they would 



obtain a fairly different set of parameters from the ones 
derived in their paper. Fortunately, the amplitudes of the 
double sine curves are not very sensitive to the adopted 
inclination angle. Considering the relative sensitivity of 
different methods on the inclination angle, the mass-loss 
rate estimation by fitting the double sine curves is more 
reliable than the one obtained by fitting the polarization 
variations near (j) = 0.5. 

3.3. Magnitudes and Distance Moduli 

For each model in Figs. |l^ and [ij, the absolute magnitude 
(My) of the system (W-R -I- O) is calculated, and placed 
in Table ||. My — —5.60 from Model B is ver y similar to 
My = -5.7 of Lundstrom & Stenholm (19841), but is far 
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from the value of Hamann fc Schwarz (199^ ), — —4.2. 
Model E gives My = —5.15 which is slightly different from 
the values estimated by Lundstrom & Stenholm (1984) 



and Hamann fc Schwarz (1992). Using the observed rriy 



4. Conclusions 

Using the 3-D Monte Carlo model developed in Paper I, 
combined with the multi-line non-LTE radiative model of 



8.27 and = 2.79 from [Lundstrom fc Stenholm (1984| ), 
we have also estimated the distance modulus (DM) for 
each model, and summarized them in the same table. The 
values of DAI from Model B and Model E are 11.1 and 10.6 
respectively. The former (Model B) agrees with the value 
given by [Lundstrom fc Stenholm (1984] ) {DM = 11.01 ± 
0.2 0), but it is slightly sm aller than the value estimated 
by [Marchenko et al. (199^ ) {DM = 11.27- 11.47). 

The absolute magnitudes of the O star corresponding 
to Models B and E are also listed in Table | They are 
M y — —5.20 and My — — 4.85 respectively. According 
to Howarth fc Prinja (1989 ), the mean absolute magni- 
tude of 06 V stars is Aly = —5.2, and that for 06 III is 
My = —6.6. Our magnitude of Model B agrees with the 
mean of their 06 V stars. Similarly, the average absolute 
magnitude of observed WN5 stars is My = —4.1 (±0.8) 



(van der Hucht 2001) which also agrees with our values 
My = -4.33 (Model B) and My = -3.58 (Model E). 

The derived O star radii are Rq = 7.2i?0(Model B) 
and i?o — 6.9i?Q(Model E). Both values are slightly 
smaller than the recent value Rq = 8. 5 ±1.0 Rq estimated 
by STL2, and are not cons istent with Rq ~ 10.0 i?Q 
( I Cherepashchuk et al. 1984[ ) and Rq ^ 8.4 - 9.3 i?0 
( Moffat fc Marchenko 1996 ) . The rather small O star radii 
of the models and the model My values compared to the 
mean observed My suggest that the O star is more likely 
to be a main-sequence star rather than a giant star. A 
similar conclusion about the O star luminosity class was 
reached by STL2. 



Vanbeveren et al. (1998) quoted Cherepashchuk 



(1975), "The observed mass of the WNE component is 
9 Mq ; its progenitor must therefore have had a mass of 
~ 30 Mq. Thus, the age of the binary is about 7 million 
years. The OB companion is an 06 star. Since the age of 
a normal 06 star is about 1-2 million years, the OB star 
must have been rejuvenated, i.e. mass transfer must have 
occurred." This argument suggests that the O star is not 
a normal O star. According to Vanbeveren et al. (1998| ), 
a quasi-conservative Roche lobe overflow model can pro- 
duce a system which resembles V444 Cyg. At the end of 
the mass transfer, when the secondary has been entirely 
mixed, the 06 component may be an over-luminous star 
with Mo — 26 Mq. Contrary to the last remark, our best 
fit models (Models B and E) suggest that the absolute 
magnitude of the O star is similar to the average of 06 V 
stars, as discussed earlier. 

Unfortunately, we can not strongly conclude that ei- 
ther Model B or E is a more favorable model based on 
the comparisons with the observed magnitudes, distance 
moduli and the results from earlier works. 



Hillier fc Miller (1998 ), we first investigated the basic be- 
havior of the WN star diagnostic lines (Hen 5412 A and 
He I 5876 A) as a function of the phase of V444 Cyg 
(Figs. ^-^). Since the emission lines from the W-R star 
envelope are contaminated by the O star atmospheric ab- 
sorption lines and the line emission from the shock-heated 
region, we first examined at which phase the profile shape 
of helium lines, especially Hei 5876 A, is least affected 
by the contaminations. We then fit the Hei 5876 A profile 
with the multi-line non-LTE radiative transfer model with 
spherical geometry to obtain a realistic model for the W-R 
star. The best orbital phase for this purpose is found to 
be around <p — 0.8. 

A summary of the derived parameters of V444 Cyg 
with our best fit models (Models B & E) is given in 
Table |[ Both models show excellent agreement with the 
observed light curves of /, Q and U Stokes parame- 
ters (Figs. Jl^ and |l^). They also fit the spectrum of 
He II 5412 Aand Hei 5876 A very well (Figs. | and |ll|). 
The effect of ignoring the presence of the bow shock in the 
continuum polarization calculation was also considered. 
This was done by first approximating the opacity in the 
bow shock region via fitting Hei 5876 A line at (f> = 0.419 
where the emission from the bow shock is most prominent 
(Fig |l3|). According to this model, the bow shock region 
has little or no effect on the continuum light curves at 
A = 5630A for aU /, Q and U (Fig. |l|). Similarly, the 
presence of the 0-star wind and the exact details of the 
bow shock are found to be unimportant for the continuum 
fits (Fig. ^6|). We also found that the spectroscopic solu- 
tion, when combined with the strong O star spectrum, is 
very insensitive to the value of the W-R core radius (c.f., 
Fig. land Fig. 

Small discrepancies are seen in the oscillation of the Q 
and U light curves (Figs. ^ and |l^) near the secondary 
eclipse ((/> = 0.5), which are possibly due to the lack of 
knowledge of the exact opacity and emissivity in the bow 
shock region. This discrepancy could be a result of the 
wrong combination of the orbital inclination angle and 
the O star radius in our model. The shape of the polariza- 
tion curves near (f> = 0.5 is very sensitive to the inclination 
angle, and is difficult to fit consistently with other obser- 
vational features. 

The most noticeable difference between our best model 
and the observations is in the IR light curves (2.2 /Ltm) in 
Figs. ^ and ^ The model correctly predicts the depth 
of the primary eclipse, but that of the secondary eclipse is 
too shallow. We argued that the difference could be caused 
by the use of the wrong monochromatic luminosity ratio 
{q) at A = 2.2 ^m. 

The mass-loss rate of the W-R component determined 
in our analysis is Mwr = 0.6 (±0.2) x 10~^ Mq yr~^, and 
this conclusion is insensitive to the model radius of the W- 
R star. The mass-loss rate of the W-R star obtained by 
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STL2, MwR = 0.6 X lO^^MQyr"^ is in good agreement 
with our value derived here. Our derived mass-loss rate lies 
between the values obtained by the orbital period change 
method (Table |l|). The fits did not allow a unique value 
for the radius of the W-R star to be derived. The range of 
the volume filling factor for the W-R star atmosphere is 
estimated to be 0.050 — 0.075 for the corresponding range 
of the W-R star radius, 5.0 — 2.5 Rq. 

Future work includes: 1. obtaining improved IR light 
curves for both models and observations, 2. fitting and 
modeling the UV spectrum and polarization, 3. checking 
the consistency between the polarization model and the 
hydrodynamical model more carefully, 4. investigating the 
cause for the asymmetry of the polarization light curve 
around the secondary eclipse, and 5. obtaining a tighter 
constraint on q. Future interferometer observations of this 
system will also prove invaluable to our further under- 
standing. 
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